# #ANALYZE DATA

##STEP 1: Install and load packages
#install.packages("tidyverse") #install tidyverse package
#install.packages("RSA") #install RSA package
library(tidyverse) #load tidyverse package
library(RSA) #load RSA package

##STEP 2: Load dataset from desktop
df <- read_csv2(file.choose()) #read "dataset.csv" from desktop

##STEP 3: Centering predictor variables
grand.M <- mean(c(df$X, df$Y)) #center x and y on grand mean across both variables
df$Y.c <- df$Y-grand.M
df$X.c <- df$X-grand.M
df$X.s <- scale(df$X)
df$Y.s <- scale(df$Y)

##STEP 4: Run RSA function 
rsa.model <-  RSA(formula = Z ~ X.c*Y.c, #Specify which variable is outcome and which are the two predictors
                  data    = df,   #Specify data
                  na.rm   = TRUE,       #Remove missing variables
                  out.rm  = TRUE,       #Default: remove outliers
                  models  = c("full"),  #Use the full polynomial equation (X + Y + x^2 + XY + Y^2)
                  missing = "listwise") #listwise deletion of missing values to be comparable to traditional methods

summary(rsa.model)  #This displays:
# 1) Percent of matches & mismatches
# 2) The polynomial model (i.e., unstandardized slopes and R^2)
# 3) Four RSA coefficients (a1 - a4)

##STEP 5: Plot Response Surface in 3D space
plot(rsa.model,           #Fitted model
     type = "3d",         #For those who have difficulty interpreting 3-d figures try "interactive" which will let you rotate the figure
     xlab = "Buyer Dependence (BD)",
     ylab = "Supplier Dependence (SD)",
     zlab = "Buyer Satisfaction (SAT)",
     surface = "predict", #Response surface is based on predicted values of outcome
     rotation = list(x = -63, y = 32, z = 15),
     legend  = TRUE,     #TRUE displays color legend
     param   = TRUE,      #Display RSA parameters
     coefs   = FALSE,      #Display polynomial coefficients 
     axesStyles = list(LOC = list(lty="dotted",  lwd=3, col= "red"), LOIC= list(lty="solid",  lwd=3, col="blue")),
     points = FALSE, #Display raw scatter points
     axes    = c("LOC", "LOIC"), #Display line of congruence and line of incongruence
     project = NULL,      #TRUE displays projections onto the bottom of the plot
     hull    = FALSE,     #TRUE displays a bag plot on the surface
     xlim    = c(-2,2),   #Set range for x axis to fit the range of value for predictor 1
     ylim    = c(-2,2)  #set range for y axis to fit the range of value for predictor 2
     )